{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "99a8eaa0-e84c-457f-8c07-36c8c5846cc2",
   "metadata": {},
   "source": [
    "# Using Tabulated Potentials\n",
    "\n",
    "OpenMM's [custom forces](https://docs.openmm.org/latest/userguide/theory/03_custom_forces.html) allow you to use potential energy expressions with arbitrary functional forms in your simulations.  However, in certain applications, you may want to run simulations with potentials defined by lists of tabulated values, rather than analytical expressions.  This can be implemented in OpenMM using custom forces' support for tabulated functions.\n",
    "\n",
    "Tabulated potentials are often used in coarse-graining; for this example, we will use coarse-grained pairwise potentials for an equimolar mixture of methanol (A) and ethanol (B) generated [[1]](#Acknowledgements) with the [OpenMSCG](http://software.rcc.uchicago.edu/mscg/) package.  You can generate tabulated values in any way convenient for your application.  Here, we load them from a text file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a95bc7a2-f6cd-4c7a-a69f-72d3c8291c6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import openmm\n",
    "import openmm.unit\n",
    "\n",
    "table = np.loadtxt(\"tabulated_potentials.dat\")\n",
    "n_entries = table.shape[0]\n",
    "u_aa, u_ab, u_bb = table.T * openmm.unit.kilocalorie_per_mole\n",
    "r_cut = 12 * openmm.unit.angstrom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bbf8d69-78fd-4071-a754-dd7ecc134d18",
   "metadata": {},
   "source": [
    "This sample file contains tabulated pairwise potentials for interacting A-A, A-B, and B-B pairs in a binary mixture.  The table values are equally spaced, with the first values corresponding to the potentials at zero separation distance, and the last values at `r_cut` (where the potentials vanish).  Note the use of OpenMM's [units](https://docs.openmm.org/latest/api-python/app.html#units) system to ensure that unit conversions are handled correctly.\n",
    "\n",
    "We can inspect these potentials by plotting them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0d1afd51-14c6-4ba5-9732-1d2f497a1220",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAGyCAYAAAAGdNXrAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAaO1JREFUeJzt3Qd4VFXaB/D/tEx674UeegcBQRQVEFTQtTfs7qoL6rquivrp2lfXsq5dRNG1YcGGDQvSpfcaICGB9N4n077nnJsJCQmQzNzJTGb+v+eZvTOTmTOXS9Z5ec97zqux2+12EBEREfkYradPgIiIiMgdGOQQERGRT2KQQ0RERD6JQQ4RERH5JAY5RERE5JMY5BAREZFPYpBDREREPkkPH2Kz2ZCbm4uwsDBoNBpPnw4RERG1g9iyr6qqCsnJydBq1cu/+FSQIwKctLQ0T58GEREROSEnJwepqalQi08FOSKD47hI4eHhnj4d6kSv/JqBN5YfxJVj0vDgeQOP/8I/XgeWPgkMuhiY+ZLTn7c+fz3m/DYHPcN74uPzP3Z6HCIiAiorK2WSwvE9rhafCnIcU1QiwGGQ41/CwsOhNQbL2wn/7qNiAKMGMJjFL4rTn9fN2g26IB0qtZX8XSMiUonapSYsPCafYNArv8oNlpO0YjOEKMeGGpc+LzowWh7LTeWw2CwujUVERO7BIId8gkGn/CqbrbYTvzAgWDk21Lr0eZHGSGg1ymeW1Ze5NBYREbkHgxzyCQE6TTuDnMZMjtm1TI5Oq5OBjlBaX+rSWERE5B4+VZND/kvflMnpnOkqISYoRgY4JXUlLo9FRL6zlUlDQ4OnT8PrGAwG6HS6Tv9cBjnkEzp7uqp5XU5JPYMcIoIMbjIzM2WgQ61FRkYiMTGxU/exY5BDPsHQ7umqUPUyOYEx8sjpKiISm9nl5eXJbIVYCq3mhna+cG1qa2tRWFgoHyclJXXaZzPIIZ8Q0N5MjiH4aE2O3S7WKzr9mczkEJGDxWKRX+Rix97g4Mb/zlCToKAgeRSBTnx8fKdNXXlVqPnPf/5TprGa30Rqi6i901UNJ6vJcRQe222AxeRyTY7Amhwislqtyn9iAgJ4MY7DEfyZzWb4bSZn0KBB+OWXX5oee6JQibruPjlmSztXVzmmrAyBTn8mp6uI6Fjsm+hd18brghy9Xs/sDTldk2M5WcGfVgfojIDV1LiMXMnGOIOZHCIi7+ZV01VCRkaGnNPs2bMnrrjiChw8ePC4rzWZTLLfRfMb+fvqqpNMVzXP5ri4wspRk8PCYyIi7+RVQc7YsWPx/vvv46effsK8efOQn5+P8ePHo6Sk7SW6Tz/9NCIiIppu7EDuv5pqck42XdUiyHFthVXz6SqxeoCIiLyLVwU506dPx8UXX4whQ4Zg8uTJ+O677+Tz7733Xpuvnzt3LioqKppuovs4+ad2LyE/doWVC6KDlEyO2WZGlbnKpbGIiDxp9erVsgZ22rRp7X7PRx99JN9z6623wlt5VZBzrJCQEBnwiCmsthiNxqaO4+w87t/avYRcxUyOUWdEqEHZd4crrIioK3vnnXcwZ84crFy5EtnZ2e1+z7333otPPvlELp/3Rl5XeHxszc3u3bsxceJET58K+WRNjusbAoq6nGpztQxyekb0dHk8IvINYgq7zqwsK+9sQQZdh1Yy1dTU4NNPP8X69etlmciCBQvw8MMPn/A9WVlZMvvzxRdfYOnSpfj8889x7bXXwtt4VZBzzz33YMaMGejWrZvcMOiJJ56QxcTXXXedp0+NusoS8o5kcsyu/8sjNigW2VXZbO1ARC2IAGfgwz955KrseuwcBAe0/+t94cKF6Nevn7xdc801MqPzf//3fycMlEQW57zzzpP1sOI98+fP98ogx6umqw4fPowrr7xSXuiLLrpIbqr0xx9/oHv37p4+NfLFmhyVmnQKxXXFLo9FROQJ8+fPl4GKIGpyqqur8euvvx739aI3l8j2ON4jVkKvWbMG+/fvh7fxqkyOmNcjcoahsU+MzQ5YbXbotJp2NOl0PciJC4qTRwY5RHTslJHIqHjqs9tr7969WLduHRYtWtS0V93ll18uMzV9+/bFwIEDm177wAMPyNuSJUvkFJdYLCTExsZi6tSp8j1PPfUUvIlXBTlErk5XObI5OrHpXyc06YwLVoKcotoil8ciIt8hpno6MmXkySyOxWJBSkpKi3oig8GAl156CVu2bGl6PjpaWVEqgpnS0tIWPbpEdmfz5s14/PHHvapTgff/DRB1YLpKaLDaEHiif8k0LSF3vSbHsVdOcT2nq4ioa7FYLHJvuueff15mYpoT27mIWp3Zs2e3eF7sW/f111/LmRfRhql5kCMWCf3www84//zz4S0Y5JBPTVe1r39VsOqZnOJaBjlE1LUsXrwYZWVluOmmm2QBcXOXXHKJzPIcG+T873//Q0xMDC699FJom/13VxDBjXiPNwU5XlV4TOQsrVYDfWMdjkUU5pxI03RVtSqrq4SiOk5XEVHXMn/+fLnx7rEBjiOTI6aqNm3a1OJ5MVX1pz/9qVWA43iPCJwKCgrgLZjJIZ/aK8dis568tYMxXDmaqlQLcsrqy2CxWaDX8v9SRNQ1fPvtt8f92ciRI9tsV7Nt27bjvkesijabzfAmzOSQ/y0jD2wMcupdb+gaZYyCVqOFHXY26iQi8jIMcsj/dj12ZHLqK1z+TLGKq6n4mHvlEBF5FQY55INBzskyOY3zzybXMznNp6wY5BAReRcGOeQzDHpN0xLy9k1XuZ7JERjkEBF5JwY55HOZHMvJpqsCI4/uk2M1q7fCihsCEhF5FQY55DMC2jtdZQw7el+F4mNmcoiIvBODHPK5TM5Jp6t0BsDQ2Inc5PqUFYMcIiLvxCCHfIbesYT8ZPvkqLyMvKl/FTcEJCLyKgxyyP8yOc1XWKlQfMxO5ERE3olBDvkMY2Mn8pPueNxi12PXMzkxQUf3yWlrh1AiIm+3evVq2T182rRp7eqw7rjp9Xp069YNd999N0wmE7wNgxzyGUa9rv1BjoqZHEdNjslqQrXZ9X5YRESd7Z133sGcOXOwcuVKZGdnn/T17777LvLy8pCZmYnXXntNNu584okn4G3YaId8LpNj6uSanCB9EEINoTLAEXU5YQHNVm8REXm5mpoafPrpp1i/fj3y8/OxYMECPPzwwyd8T2RkJBITE+X9tLQ0zJw5s1UzT2/ATA75YJBj7dRMTvNsTkldiSrjEVEXJ6auG2o8c+vgtPnChQvRr18/ebvmmmtklqYjU+/79u3D0qVLMXbsWHgbZnLIZwR4qCbHscIqqzILhbWFqoxHRF2c2Gz0qWTPfPYDuUBA4zYZ7TB//nwZ3AiiJqe6uhq//vorJk+efNz3XHnllbKGx2KxyFqc888/H3PnzoW3YSaH/HS6St1MTnxwvDwyyCGirmTv3r1Yt24drrjiCvlYFBJffvnlskZH1OaEhoY23Z566qmm97344ovYsmULtm7disWLF8tszqxZs+BtmMkh/8zkqNy/KiE4QR4LagtUGY+IujhDsJJR8dRndyCLY7FYkJKS0vScmKoyGAx46aWXZCDjEB0d3XRf1OP06dNH3hfTXFVVVTK7I4qPHc97AwY55HOrq9qXyYl0T5BTwyCHiOQ66w5NGXmCxWLB+++/j+effx5Tp05t8bOLL75Y1urMnj27XWOJqSuhrq4O3oRBDvnndJXKNTmOIIfTVUTUVSxevBhlZWW46aabEBHROIXf6JJLLpFZnuMFOeXl5XIlls1mQ0ZGBh577DH07dsXAwYMgDdhTQ753HSVJ1ZXJYQoQU5+bb4q4xERudv8+fNlcfGxAY4jkyOmqo63LPyGG25AUlISUlNT5TTVoEGD8MMPP8iaHm/iXWdD1AX3yWmeyRG7HltsFui1/L8WEXm3b7/99rg/Gzly5HGXkXelnd2ZySGfEeDsjscq/B82OjAaOo0ONruNe+UQEXkJBjnk3zU5dquyn4WLdFpdUzdyrrAiIvIODHLIB5eQt6MmR6x60OjcMmXFIIeIyDswyCH/zOSI5Z0q75Xj2BCQy8iJiLwDgxzyz80Am9flcBk5EZFPYpBD/rkZoHyDupmcxBClIy+XkRMReQcGOeQzjAYnMznc9ZiIyCcxyCGfEaDrwGaAApt0EhH5NAY55DMCGzM5HZ6uUqsmJ+To6qqutFkWEZGvYpBDPiOgsUGcp6ar4oOU1VVmmxllpjJVxiQiIucxyCGfq8lpdyZH5dYOBp1B7nwssFEnEXUF119/PTQaTdMtJiYG06ZNw7Zt21R9j6cwyCGfq8mx2uywWDvY2kElTRsC1hSoNiYRkTtNmzYNeXl58vbrr7/KJpvnn3++6u/xBAY55HOZHKGhPUGOykvIm9fl5NXkqTYmEZE7GY1GJCYmytvw4cNx3333IScnB0VFRaq+xxPYKpl8LpPjqMsJDjjJG4KVqSXUlap2DimhKfKYW5Or2phE1PWIxQd1ljqPfHaQPkhOIzmjuroaH374Ifr06SOnodz1ns7CIId8hl6nhU6rkdNV7arLCVEaaqKmWLVzSA5JlscjVUdUG5OIuh4R4Iz9aKxHPnvtVWsRbAhu9+sXL16M0NBQeb+mpgZJSUnyOa1Wq+p7PMG7zoZIrf5V5nYEOcExqgc5KWGNmZxqZnKIqGs488wzsWXLFnlbu3Ytpk6diunTp+PQoUPyKIIZcRs0aFC73uNNmMkhn+tfVdtgRYPV2v5MjrkGaKgFAtr/L5+TTVcdqWYmh8ifiSkjkVHx1Gd3REhIiJxqchg1ahQiIiIwb948vP3226irU6bdDAZDu97zxBNPwFswyCGfzOTUtyeTYwwDdAGAtQGoLQYCurn8+cmhynSV2Cen1lzboZQxEfkOURPTVf//r9Fo5LSTCG5SUlI6/B5vwiCHfLMTeXtWV4nCPJHNqTwC1BQBka4HOeEB4QgLCENVQ5XM5qRHpbs8JhGRO5lMJuTn58v7ZWVleOWVV2Qx8YwZM1R9jycwyCHf7ETenkyOEBLbGOSUqHYOqaGp2F26W9blMMghIm/3448/ysJhISwsDP3798dnn32GSZMmqfoeT2CQQz45XdWuTI4QHKscRSZHJWLKSgQ5h6sPqzYmEZE7LFiwQN7c/R5P4eoq8snpKpO5nZ3IHcXHoiZH7b1yuMKKiMijGOSQby4ht3RgusoNmRyBK6yIiDyLQQ75lAB9BzuRNwU5xarW5AjM5BAReRaDHPLvTE5TTU6x6pkc1uQQEXkWgxzyzSXklg7W5Kg4XeWoyRHLyCsbKlUbl4i6Rs8q8p5rwyCH/Lwmx1F4rN4ScrEBWJQxSt7nlBWRf9DpGqfKGxo8fSpeq7a2ttXOye7GJeTkm0vI2x3kOPpXFYl/ZigbBKqUzRG7Hovi4/7R/VUZk4i8l16vR3BwMIqKiuSXuLc1qvR0BkcEOIWFhYiMjGwKCDsDgxzyzc0AO5rJsdQDDTWAUemqq0Zdzo6SHThcxb1yiPyBaGsgNsfLzMz0uiaV3kIEOImJiZ36mQxyyEenq9pZkxMQAohmdpY6JZujUpDTLVxpEZFTlaPKeETk/QICApCens4pqzaI7FZnZnAcGOSQjxYetzOT48jmVGQrK6yie6pyHt3Du8tjVmWWKuMRUdcgpqkCAwM9fRrUiJOG5FM6XHjcvC5HxV2Pe4T3kMesCgY5RESewiCHfIrTmRyVl5E7gpyC2gLUmpUVBURE1Lm8Nsh5+umnZSHXXXfd5elTIV8uPHbThoCRgZGIMEbI+6zLISLyDK8MctavX4+33noLQ4cO9fSpUFdt0NmhTI76QY7AuhwiIs/yuiCnuroaV199NebNm4eoKGVDteMxmUyorKxscSP/1uHVVW7qRN58yupQJZeTEhF5gtcFOX/9619x3nnnYfLkye2a0oqIiGi6paWldco5kvdPV3WsJkf9TuQCgxwiIs/yqiDnk08+waZNm2Tw0h5z585FRUVF0y0nh3uS+DvnpqvilWNVgXumq7jCiojII7xmnxwRoNx5551YsmRJu/cYMBqN8kbU9DvhTJAToTTUROVhtwQ5mZWZcltzUUhPRER+mMnZuHGj7GsxatQo2QNE3JYtW4b//ve/8r7V2oEaC/JbHe5CLkSkKsf6CqBevboux67Hoht5ualctXGJiKiLZXLOPvtsbN++vcVzN9xwA/r374/77rvPI9tBk59kcoxhQGAkUF8OVBwGAgeqci5B+iAkhiQivyZfFh9HBZ64kJ6IiHw0yAkLC8PgwYNbPBcSEoKYmJhWzxOpuhmgEJF2NMhJUCfIcRQfiyBHtHcYHj9ctXGJiKgLTVcReWwzQCGycWVeRY576nIqMlUdl4iIulAmpy2///67p0+Buuh0VcczOaluCXJ6R/aWxwPlB1Qdl4iITo6ZHPLZzQDFiqYOTVcJYrpKRemR6fKYUZah6rhERHRyDHLIJ6erbHbAIv6nw5kcdYOcPpF95DG3JhfVDdWqjk1ERCfGIId8svC4w1NWjkxOubrTVaJRZ1yQ0jbiQAWnrIiIOhODHPLZIKdDxceOwuOqXMBqUfWc0qM4ZUVE5AkMcsin6LQa6LUaJ/pXxQNaA2C3AVV5bpmy2l++X9VxiYjoxBjkkM9xqhO5Vnu0vYPKK6wcQQ6Lj4mIOheDHPI5Lm0I6Ibi475RfeWRmRwios7FIId8jtMbAjYFOepmcnpG9IQGGpTWl6KkrkTVsYmI6PgY5JDPMRqcmK5qXnys8gqrYEMwUsOUJeoZ5dwvh4ioszDIIZ8ToHOiSacb98ppving/jIWHxMRdRYGOeTDmRwng5zybNXPqU+UUny8r2yf6mMTEVHbGOSQz2ZyOlx4HN1LOZZlqr5XzsBopbP5rpJdqo5LRETHxyCHfI7zhcfdAH0QYG0Ayg+pek6DYgc1rbCqt9SrOjYREbWNQQ75HKeXkIu9cmKV2hkU7VH1nBKCExAdGA2r3Yq9ZXtVHZuIiNrGIId8TmBjTU69uYOrq4S4fsqxSN1ARKPRYFCMks3ZWbxT1bGJiKhtDHLI5wQZdM4HObGNQU6x+gXCjikr1uUQEXUOBjnkc4IClCCnrsF7MjlCUyanhJkcIqLOwCCHfE5gYyanzpXpKpHJsdtVPa+BMcoKq4MVB1FrrlV1bCIiao1BDvnsdJVTQY5YRq7VAw3VQOURVc8rPjgecUFxsNltLD4mIuoEDHLI57hUk6MzHN0vx51TViw+JiJyOwY55HNcqslxc13OwFhlympHyQ7VxyYiopYY5JDPcakmp8UKK/WDnGFxw+RxS+EW1ccmIqKWGOSQzwl2ZHLMHdwMsFUmZ59bghytRosj1UeQX5Ov+vhERHQUgxzy3cLjBouLQc5u1VdYhRhC0D+6v7y/uXCzqmMTEVFLDHLI5wQ2ZXKcrcnpD2gNQF2ZWzqSj4wfKY8bCzaqPjYRER3FIId8OJPjZJCjNwIJyioo5G6C2kYmKEEOMzlERO7FIId8eAm5kzU5QvII5Zir/pTSiHhl7IyyDFQ2VKo+PhERKRjkkO8uIXd2ukpIUbItOKJ+Jic2KBbdw7vDDjtXWRERuRGDHPI5Lk9XNc/k5G0FbC5khE5Sl7OpQP0gioiIFAxyyKf3ybE7uzoqbgCgDwRMlUDpAbfV5Wwo2KD62EREpGCQQz47XSWYLE5mYXR6IHGo2+pyxiSOkccdxTtQY65RfXwiImKQQz48XaXalJUbgpzk0GSkhaXBardyKTkRkZswk0M+R6fVIECv/GrXemnxcfNsztq8tW4Zn4jI3zHIIZ+kevGx1Qy1jUsaJ48McoiI3INBDvn4XjkuBDkx6UBgJGCpA/K3QW2nJJ4ij3vL9qK0vlT18YmI/B2DHPJJquyVo9UC3ZRsCw6tgdpigmKQHpUu76/PX6/6+ERE/o5BDvn2MnJXpquEbqcqx2z1gxxhbOJYeeSUFRGR+hjkkE8KMmhdz+QcG+So3JFcGJukBDnr8tepPjYRkb9jkEM+PV3lUk2Oo/hYbApYWwIUZ0BtoxNGQ6fR4VDlIeTX5Ks+PhGRP2OQQz5JldVVgj4ASBmt3M9eDbWFBoRiUIzS8fyPvD9UH5+IyJ8xyCGfb+3gsu6nuq34uMWUVR6nrIiI1MQgh3xScON0Va2rmRzBscIq271Bjig+drrXFhERtcIgh3ySKvvkOKSNBTRaoPwQUJ4DtQ2LG4YAbQAK6wqRWZmp+vhERP6KQQ75pEDHPjlqZHKMYUByY4uHzOVQW6A+ECPild2VOWVFRKQeBjnk24XHamRyhJ6nuy3IEcYksY8VEZHaGOSQT1I9yOl1hnLMXOb2/XJsdpvq4xMR+SMGOeSTVNsnp3ldjs4IVOW5Zb8csYw8xBCCyoZK7Cndo/r4RET+iEEO+STV2jo4GIKAtDFHszkq02v1cmNAgS0eiIjUwSCHfJLq01XHTlm5cyl5/lq3jE9E5G8Y5JCPBzkq1rf0nKQcM1cANhWDp0ZjEpVM0aaCTTBbzaqPT0TkbxjkkG/X5Kg1XeXoYxUQBtSXA/nbobb0qHREB0ajzlKH7cXqj09E5G8Y5JBPBzm1Zot6g+r0QI8Jbpuy0mq0OCXxFHmfdTlERK5jkEM+3qBT5eXYPRvrcg66ty6HzTqJiFzHIId8kqptHdraFFD0sbI0qDs2gHGJSp+sbcXbUGuuVX18IiJ/wiCHfHq6SqyuUrXpZfxAIDgWEAHIkQ1QW2pYKpJCkmCxWbC5cLPq4xMR+RMGOeTT++RYbXaYrSoGOVrt0WyOG6asNBoNl5ITEXlDkGM2m5GTk4O9e/eitLRUrXMiUm26SvW9clrsl+OmPlaNS8lZfExE1MlBTnV1Nd58801MmjQJERER6NGjBwYOHIi4uDh0794dt9xyC9avX+/Uybz++usYOnQowsPD5e3UU0/FDz/84NRY5N8MOg10Wo1763IOrwcaatxWfLy7ZDcqTBWqj09E5C86FOS8+OKLMqiZN28ezjrrLCxatAhbtmyRmZw1a9bgkUcegcViwZQpUzBt2jRkZHSsx09qair+9a9/YcOGDfImPuOCCy7Azp07O/rnIj8npn2OrrBSOciJ6glEdANsZiBrlbpji7Kf4Hj0iugFO+zYkK9+3Q8Rkb/Qd+TFq1evxtKlSzFkyJA2fz5mzBjceOONMiPzzjvvYNmyZUhPT2/3+DNmzGjx+Mknn5Rj/fHHHxg0aFBHTpVI1uVUmyzqT1dpNEDvM4FN7wEHfgP6TnXLlNXBioNyKfnZ3c9WfXwiIn/QoSDns88+a9frAgMDcfvtt8MVVqtVfl5NTY2ctmqLyWSSN4fKykqXPpN8S7BjQ0C1MzlCn7Mbg5xf1R9bLCVPGodP9n7CPlZERJ0V5Nx9993tfu0LL7zgzPlg+/btMqipr69HaGgovvzyS1nz05ann34ajz76qFOfQ77PbXvlOOpyNFqgeB9QngNEpqk6/OjE0dBAg8yKTBTWFsopLCIicmOQs3nz5nbXQzirX79+ss6nvLwcX3zxBa677jo57dVWoDN37twWgZfI5KSlqftlQ11XoGOvHHdkcoKigJTRwOF1ypTVqOtUHT7CGIEBMQOwq2SXXGU1o3fLqVwiIlI5yBH1OO4WEBCAPn36yPujR4+WK7VeeukluaLrWEajUd6I2hJkUOrqVa/JaT5l5aYgx7HKSgQ56/LXMcghIvLFzQDFbrXN626I2qtpdZW7gpzeZynHg78DNvU/Y2yispRcZHJU3bWZiMhPdCiTcywxpTR//nzs3r1bTlENGDAAN910k9w/xxkPPPAApk+fLqecqqqq8Mknn+D333/Hjz/+6Mppkp+3dnBLTY6QPBIIjADqy4Ejm4A0pYO4WkbEj4Beq0deTR5yqnLQLbybquMTEfk6pzM5Yh+b3r17y71zxG7HxcXF8r54btOmTU6NWVBQgFmzZsm6nLPPPhtr166VAY7Yd4fI2dYOblldJej0R7uSiykrlQUbgjE0dqi8vzZ/rerjExH5OqczOX/7298wc+ZMuTGgXq8MIzYCvPnmm3HXXXdh+fKOb3kvskJEagkzKr+XNSaL+y6qqMvZ/Y2ylHzSfW5ZSr6pcJOcsrq076Wqj09E5MtcyuTcd999TQGOIO7fe++98mdEnhYaqPxuVtW7Mchx1OUc3gDUlbutxcO6vHWw2W2qj09E5MucDnJEb6ns7OxWz4uGnWFhYa6eF5HLQo0GeRS7HrtNZDcgJh2wW93SsHNI7BAE6gJRZiqTe+YQEVEnBDmXX365LDJeuHChDGwOHz4sC4XFdNWVV17p7LBEqmdyqt2ZyXFMWbmpLsegM2BonFKXI6atiIioE2pynnvuObmi6tprr5W1OGKJq9jj5rbbbpNNNok8LdSoc38mxzFltfYNpS5HLPV2YTPM462yEnvlbC7YzLocIqLOCHJEQCM26ROtFQ4cOCCDHLGJX3BwsLNDErlluqrK3UFOj9MAXQBQng2UHABilc0s1TIyfqQ8MpNDRNSJ++SI/lI7duxAYWEhbDYbsrKymn4mVl4ReVJo4+qq6nqzez8oIAToNk6pydn/i+pBjpiu0mq0OFJ9BAU1BUgISVB1fCIiX+V0kCP2rxF72pSUlLT6mZjGEl3EiTwpzFGT4+5MjtBnSmOQ8zMw7lZVhw4NCEW/qH7YXbobm4s2Y1rINFXHJyLyVU4XHs+ePRuXXXYZ8vLyZBan+Y0BDnlXJqcTgpz0qcoxcwXQUKv68KIuRxB1OURE5OYgR0xRiQ7gCQlMnZN3r66qabDCanNz76e4fkBEGmA1AVkrVR9+REJjkFPIIIeIyO1BziWXXCL7ShF5eyZHqGlwczZHrKhKb2w/krFE9eFHxClBzt6yvahuqFZ9fCIiX+R0Tc4rr7yCSy+9FCtWrMCQIUNgMCgrWRzuuOMONc6PyGlGvRYGnQZmq11OWYUHtvwddcuU1YZ3lCBH5aXkotg4JTRFFh9vK9qG8SnjVRubiMhXOR3kfPTRR/jpp58QFBQkMzqi2NhB3GeQQ54mfg9FNqes1tw5xcc9T29cSn4IKNkPxKarXpcjghyxlJxBDhGRG6erHnroITz22GOoqKiQS8czMzObbgcPHnR2WKKu17+q+VLy7hPcN2XlKD5mXQ4RkXuDnIaGBtnaQat1eggi3+hf1Zwb63IcmwKK6Sqzzc17/xAR+QCnI5TrrrtO9q0i8mZhnbmMvPlS8kOrAZO6BcK9InshPCAc9dZ67CnZo+rYRES+yOmaHLEXzrPPPivrcoYOHdqq8PiFF15Q4/yI1GnSaeqkzEdMHyCyu1KXIzYH7H+uakOLXY/FlNWyw8tkXc6QuCGqjU1E5IuczuRs374dI0aMkNNVorXD5s2bm25btmxR9yyJXFxG3ik1OU1LyRuzOWL3Y5WxLoeIqBMyOUuXLnX2rUQeyOR0UpAjiCBn/Twg42fVl5KPTBjZVHwsmuI2X9VIREQuZnIeeOABrFu3rqNvI/LJmpwXluzFOS8ux9I9hcd0JTcCFTlAkbq1M4NiBiFAG4DS+lJkVR5tiEtERCoEOaJX1fnnn4+kpCT8+c9/xnfffQeTydTRYYg6t3+VGzI5RVUmvPb7AewtqMINC9bjsW93wSbaRwQEAz0nKi8S2RwVBegCMDh2sLy/pZDTwkREqgY57777LgoKCvDpp58iMjISf//73xEbG4uLLroICxYsQHFxMa84ed8+OW4Icr7YdBgWmx0RQUrR/TurMrFkV/7RruRu3i9HFB8TEZHKhceiDmDixIlyddWePXvk9NW4ceMwb948pKSk4PTTT8dzzz2HI0eOODM8kdd3Ihf1MAvX58j7D5zbH7PGdZf3l2cUt9wvJ3sNUF+p6mc3r8shIqLjU2UnvwEDBuDee+/FqlWrkJOTI/fQET2tPv74YzWGJ3JamJsKj9dmliKzuAYhATqcPzQZE9Nj5fN/HCxRXhDTG4juDdgsQOYyVT97WNwweTxUeQjFdcycEhEdj+rbFcfHx+Omm27C119/jXvuuUft4Ymc2/FY5UyOI4szc3gyQox6jO0ZIxdRHSyqQUFlvVt3P44wRqBPZB95n3U5REQqLiEXtTcno9frkZiYiClTpmDGjBkd/Qgir15CLqaqftldIO9fMipVHiOCDRiUHI4dRyplNueC4SlKkLP2DSDjF9WXko9KGIX95ftlXc7k7pNVG5eIyK8zORERESe9ic7kGRkZsrfVww8/7J4zJ/LQ6qrs0lq5uWCATouhqZFNz4/rGdNyyqr7aYA+CKjKBQp2umdTwALW5RARqZbJEaur2kssL7/ttttkt3IiT9fkqLV5nsjWCP2TwmDQHf13wqm9Y/D2ykz8cbBUecIQCPQ8Hcj4SbklKku/1WzWubt0N2rNtQg2BKs2NhER/L0m54MPPjjuz/7xj3/I44QJEzB69GhnP4JItUyO1WZHvdmmyhXdkVshj4NTIlo8f0rPaGg1kAXJ+RWNdTl9z1GOe76DmpJCk5AYkgir3YptxdtUHZuICP4e5MyePRuLFy9u9fzf/va3pgBI7KOzaNEi186QyAXBAbqmUpgqlZp07jjSGOQktwxywgMNTYHP2szGKav+54lNF4AjG4EKdbdU4JQVEZGbgpxPPvkE11xzDZYvX9703Jw5c+QmgexrRd5CTE+FBqi3V46Y8tqZq0xXDU4Jb/VzR5CTUVCtPBGWCKSNdUs2xzFlxU0BiYhUDnKmTZuGN954AxdeeCE2bNiA22+/XWZtRIDTv39/Z4cl8uoVVnkV9SitaYBeq0HfhLBWP+8VGyKPB4sbgxxhQOMKw93fwB2ZnK1FW2ER+/EQEZE6XciFK664AmVlZTjttNMQFxeHZcuWoU8fZf8OIl/c9dgxVZWeEIZAg67Vz3vHhcqj2C+nyYDzgSUPAodWATUlQIiyCstVYq+cMEMYqsxV2Fu2VzbvJCIiJ4Ocu++++7gbAI4YMQKvvfZa03MvvPBCR4Ym6hL9q3Y4pqqSW09VCb3ilEyOKD4WzTq1ohI5qgeQOATI3w7s/R4YOQtq0Gl1GBY/DCuPrJRLyRnkEBG5EORs3tz2nhy9e/dGZWVl08/VWKZL5I2ZnJ1H2l5Z5ZAaFQyDTgOTxYYj5XVIi25c2j1gphLk7FmsWpDj2BRQBDmiLueagdeoNi4Rkd8FOSwoJn/vX3V0+XjbmRydVoPuMSHYX1iNg8U1R4Oc/ucDS58EDvwGmKoAY+t6HpdWWBVuVm0fICIiX6F67yoiX931uLLejIJKk7zfVtFxq+LjombFx/EDlIad1gZVe1kNjh0Mg9YgG3Uerjqs2rhERH4X5GRnZ3do8CNH1N0XhMgZYYFKk87KOtf2yclsLCaOCzM2jdmWXm0VH4sMS9Mqq9b7SznLqDM21eJwKTkRkQtBzimnnIJbbrkF69atO+5rKioqMG/ePAwePJgbAZJXiApWAhKx9NsVophY6NmYqTkeR/Fxi2XkjrocQWRyzI07IqtgRIIyZcUgh4jIhZqc3bt346mnnpJ75BgMBtmyITk5GYGBgXIp+a5du7Bz5075/L///W9Mnz69I8MTuUV0iFEey2pdC3JEjY3QuzGIOR7Hzx2ZnybJI4CwZKVh58HfgX7ToNamgO/iXWwq2KTKeEREfpnJiY6OxnPPPYfc3Fy8/vrr6Nu3L4qLi2XHceHqq6/Gxo0bsWrVKgY45DWiQ9TJ5DhqbE6ayYlVpqtyK+pR29CsDkirVfbMUXljwOFxw+UxqzILpfWNzUGJiMi5zQBF5uaiiy6SNyJvFxUcoPJ0VeiJPy8kAJHBBpTXmuV7BjXvcSWmrNa9pSwlt/wH0Cvn5orIwEj0juiNAxUH5Cqrs7ud7fKYRES+gKuryOfFhLoe5Ijl2e2tyWm5wuqYKavu44HQBKC+QpmyUrkuR2wKSERECgY55DeZnMp6C8xWm1NjFFaZUNtglfvgdHPsfXMCjhVWWY2BUROtDhh4gXJ/55dQu1mnyOQQEZGCQQ75vMjgALmC25XiY0dGJi0qCAH6k//fJjUqSB7FrsetDPrT0a7kFmXfHVeNTFCCnF0lu1BnaeMziYj8EIMc8nki+xIZpBQfl9U4t1eOYzl4e6aqHO0dhMNlbQQcaeOAsCTAVAEcWAo1JIckIz44Hha7BduLtqsyJhFRV8cgh/yCKAZ2pS7HsRzcMQ11MimRJ8jkiFVWAy9UdcpKtHNwTFltLNyoyphERH4d5JjNZuTk5GDv3r0oLeXSVfJeMa4GOR0oOm4xXVVWJ7uRn3DKSqWNAUcnjJbHDfkbVBmPiMjvgpzq6mq8+eabmDRpEiIiItCjRw8MHDgQcXFx6N69u9wRef369e45WyJXl5HXuhbkOFZNnUxiRCC0GqDBakNxdRt1N6mnAOEpQEMVcOBXqOGUxFPkcWvRVjSIHllERH6uQ0HOiy++KIOat99+G2eddZZs27BlyxaZyVmzZg0eeeQRWCwWTJkyRe6K7NgkkMhrlpFXd/zLX6zIyi6tlfd7nmS3YweDToukCCWbk1PWOVNWPSN6IjowGiarCduLWZdDRNShIGf16tVYunQp/vjjDwwaNAgTJkzAkCFD0KdPH4wZMwY33ngj3n33XRQUFGDmzJlYtmwZrzB5VSbHmdVVOaW1sNjsCDLokBge2O73nbAup/mU1d4fAHOdKnU5jimr9fnMphIRdSjI+eyzz2RQo9PpcM0116CoqKjN1xmNRtx+++24+eabeYXJK0S7UJPTvB5HBBLt5ajLOVxWe5wXjAYi0oCGamD/L1BzympDAetyiIicLjwWmZvMzExeQfKfIKedU1UOKc2Kj9skAqZBjVNWOxZB1bqcQtblEBE5HeTccccdeOCBB+TqKiJfXkJ+oKhjRcetMzknmIoa2DhllbFElY0Be0X0knU59dZ67Cje4fJ4RER+GeRceumlchWVqM0RU1eiGFl0IG9o4KoO8j7RLtTkZDZuBNiro5mcyOAT1+QIySOA0ERlyipzBVwlptNGJYyS9zllRUT+zukgR0xVffnll7jnnntQW1uLp59+Wk5hhYaGYujQoeqeJZFK01UlNQ2y2aY7uo+fqCbnuJ8pVln1P1e5v/c7qDllxeJjIvJ3emffKPbEEbcLLmhsNgigqqpKLinftm2bWudHpGqQ02CxyUabIcb2/erXmCwoqDR1aCNAh6RIZSVWvdkmp8liQo1tv7DfecCGd4A93wPnPq8EPi5wrLAS++WYrWYYdEpLCyIif6NqW4ewsDBMnDgRf/3rX9UclshlwQE6GBsba3akLseRxYkNDUBEY/+r9jLqdUgIN568LqfnRCAgDKjOB3Jd7yLeO7I3ooxRslHnzpKdLo9HROQXQU52dnaHBj9y5EhHz4fILUStijMrrA52sJ1Dh/fKEfRGoM/Zqk1ZaTXaprocTlkRkT/rUJBzyimnyLYN69atO+5rKioqMG/ePAwePFjuiNwRoq5HfIbICMXHx+PCCy+UuykTeaq1g6Mxp7NBjqMb+XGXkTv0P085iikrFYxObOxjxf1yiMiPdagmZ/fu3XjqqadkywaDwYDRo0cjOTkZgYGBKCsrw65du7Bz5075/L///W9Mnz69QycjdkgWU10i0BHtIR588EFMnTpVjhsS4tyXDJErrR0cK6s6WnR87F45x90Q0KHPZECjBYp2A+U5QGSaKsXHmws3w2wzw6BlXQ4R+Z8OBTnR0dF47rnn8MQTT+D777/HihUrkJWVhbq6OsTGxuLqq6/GOeecI7M4zvjxxx9bPBYtIkRGRyxNP/3001u93mQyyZtDZWWlU59L/sGZ1g5NjTk7uHz82Omq3IqTdBoPjlaaduasBfb/DIy+Ea7oE9kHEcYIVJgqsKtkF4bFDXNpPCIiv1ldJTI3ogv5RRddBHcSU1+O4Op401uPPvqoW8+BfEdHa3LEsu+DTm4E2Kom52TTVUKfKUqQk/GLy0GOqMsRq6x+zf5V1uUwyCEif+T06iqRuenWrRtmzJiBhx56SPa12rdvX4f3IDkeMc7dd9+N00477biZoblz58pAyHHj7st0IjGNQU5RVft2Fi6qNqHKZIFWA3SLUWprOiq5KZPTjiAnfYpyzFwGWFzfVJP75RCRv3N6nxxRJyP2xNm8ebPc+fjNN99EaWkpgoKC5C7Ia9eudenEZs+eLffbWbly5XFfIxqBihtRx+pj2tfxe3+BUo/TPSZELgd3RnLjXjnltWa5584J9+dJHAqExAM1hUD2GqDXGXDF2MSx8ripYBP3yyEiv+R0Jqd///644oor8Mwzz+Cnn35CYWEhFi9ejMTERJx9duNyWCfNmTMH33zzDZYuXYrU1FSXxiJySItWsjE5JysCbrSvoEoe+8Q7V3QshAUaEB6oBDZ5J8vmiE0ARQGyIOpyVNgvx9HHalsxN+gkIv+jVXMfErGa6oMPPkBubq7TU1QigyOWnv/222/o2bOnWqdHhLTG5dx5FfWwWG0nvSL7CpVMTt8E54Oc5lNW7cogpTcGOaIuR4X/T45JHCPvr8s7/rYPRES+yukgx2Zr+0ti3Lhx+P33350aUywfF0HSRx99JPfKyc/PlzexeovIVfFhRgTotLDa7DLQae90VXp8mEuf27TCqvzkn4leZ4rwRFlKXpkHV41JUoKctfmuTR8TEflVTY5oxCkKgocPH45hw4bJY79+/eRGgdXVypdDR73++uvyKFZuHbuU/Prrr3f2VIkkrVYj63LEsnAxZeWYvjpeVnFfoTJdla5SJif3RLseN19Knjxcae8gCpCHXaFKXY7oYyXaPATplXMhIvIHTgc5Ykpp69at8vbqq68iIyNDZndEivzxxx93aky1VmYRnagzuAhyDpfWAb1PvLJKFAuLlVW940JVKXhuV5Aj9JqkBDkHf3c5yEkLS0NiSCLya/LlxoDjk8e7NB4RkV8EOWLXY3FzqK+vx4EDBxATEyOLj4m8kSN7c7IdiB1TVd2igxFocG5lVauanI4EOStfBA4uE5G/KK5xuS7nmwPfyLocBjlE5E9UKzwWGwSKpeMMcKgrFB/nnKQI2LGyKj3BtXocIaVxGXm7Mzlp4wB9IFCVCxRnuPz5Y5OUKSs26yQif6NakEPUVaarhJzS2natrEp3Yfn4sZmc/Ip6WfR8UoZAIE0JTOSUlYscK6x2lOxAVYMSvBER+QMGOeRX2rtXTkZjJqevCpmc+LBA6LUaWGz2du+2LKesVApyRE1Ot7BusNltcmNAIiJ/wSCH/EpaYyanoNIEk8V6/JVVjuXjLq6sEnRaDRIjlCmrI+W1HQtyslYCtrbPsyO4lJyI/BGDHPK7Jp3BAboTNs0UK6sq6tRZWXXslNWR9uyVIyQNA4wRgKkCyHd9t2LHUnJuCkhE/oRBDvkVsdqoqS7nOEHOluxyeewVF+ryyqrWGwK2s/hYqwO6jVPuH1rt8uePThwtj3vL9qKsvszl8YiIugIGOeS/K6yOU3y88ZASBIzuHqXaZzoadbY7yBF6TFCOWatc/vzYoFj0iewj73OVFRH5CwY55Md75bQdcGxoDHJGqRjkpEQGn3CKrE3dT1OO2atFHxXVlpKvy2cfKyLyDwxyyO84pqsOFLVuP1JvtmL74Qp5f3SPaNUzOUc6kskRdTkBoUBdGVC4y+VzaGrWySCHiPwEgxzyOyO6KRma9VmlsB2zb82OIxVosNoQExKAHjHH723l9pocQac/ul/OoVWq1OVoNVpkVmSisLbQ5fGIiLwdgxzyO0NTI+QKK9Gbak9+1XGnqkSRslocq6sq6y2oqjc7UZez0uVzCA8IR//o/vI+szlE5A8Y5JDfMei0OKVxKmrNwZK2i457qFePI4QY9YgMNsj7ue1dRi50n3B0hZUKDWy5lJyI/AmDHPJLp/aOkcc1B4pbbAK4qSmTo149jkNyhBNTVskjlT5WtcVAyX7VNgVkJoeI/AGDHPJLp/ZSgpy1maVN/aQyCqtRUtOAAL0Wg1PCVf/MoxsCdiDI0QcAKaOU+9l/uHwOI+NHQq/R40j1ERyuOuzyeERE3oxBDvmlQcnhCDPqUVVvwc5cZTXV2ysOyuPp6bEw6tXZBNClbuQOjuLjHNeDnGBDMIbEDZH3mc0hIl/HIIf8kl6nxZiejXU5B0rkxoCLNh2Rj28/U9k0T20pUU5kcgTHzscqZHIELiUnIn/BIIfg73U5by4/iHs/3ya7hE9Mj8XIxiXm7pqu6nAmJ/UU5ShqcmqO1hC5vClg3jpZh0RE5KsY5JDfunR0GvonhqG0pqFpldWcs9Ld9nlHg5wOrK4SgqOBOGXpN3LWunweQ+OGwqgzoqiuCJmVmS6PR0TkrRjkkN+KCDLgq79OwE2n9ZSPz+4f3zSF5Q6ODQHzK+thsdqcq8tRYcpKBDjD44bL++xKTkS+jEEO+TXRZfz/zh+IjQ9NxuvXNK5icpO4UCMMOo1czVVQZXKuLkeFTI7ApeRE5A8Y5BABiAk1yqXjbv0/m1aDJGf2ymmeycndDJg7ON11kuJjm9315p9ERN6IQQ5RJ3I06uxwkBPdCwiJA6wNQN4Wl89jUOwgBOuDUWGqwL6yfS6PR0TkjRjkEHUipzYEFEQfLRXrcgxaA0YlKNNza/PUmQIjIvI2DHKIOlGqI8gp62CQ44a6nKal5PnrVBmPiMjbMMgh8kAm57AzQU5asyBHjWadjUHO+vz1aBDTYEREPoZBDlEnSosOlsecstqOvzlpWGOzzhJVmnX2i+qH2KBY1FnqsLlws8vjERF5GwY5RJ0oLSq4KZNja2wM2qFmnaIruUp1ORqNBuOTx8v7q46scnk8IiJvwyCHqBMlRQZCqwEaLDYUVXdwr5wWdTnq9LE6LeU0eVyZu1KV8YiIvAmDHKJOZNBpm/bKEU1BO6ypWac6xcenJp0KDTTIKMtAYW2hKmMSEXkLBjlEnSwtOsj5upymZp0ZQI3Sb8sVkYGRGBw7WN7nlBUR+RoGOUSdrJuj+LjUiRVWKjfrFCakTJDHVbmsyyEi38Igh8hDxcdOTVfJAcaqWpczIVkJctbkroHFZlFlTCIib8Agh8hDy8iznQ1yVK7LEdNV4QHhqGyoxI7iHaqMSUTkDRjkEHmoJsepDQFbNOvcpEqzTr1Wj1OTT5X3OWVFRL6EQQ6Rh6ar8irqYLY60QFc5WadzaesWHxMRL6EQQ5RJ4sLM8Ko10LsBdjhbuRuaNYpODYFFNNVZfVlqoxJRORpDHKIOpnYaTg1yrFXjpNTVio360wISUB6VDrssOOPPHUCJyIiT2OQQ9TVeli5oVmncFpy4+7HR7j7MRH5BgY5RF1xGbnKzTpb7JdzZBVsdidqhYiIvAyDHCIPrrByehm5ys06hZHxIxFqCEVJfQm2F29XZUwiIk9ikEPkAd2iQ1wLcuQg6m4KaNAZmhp2/p7zuypjEhF5EoMcIg/oEatMV2UW18DubE1NmrqbAgqT0ibJ49LspaqNSUTkKQxyiDyge2Mmp6regtKaBucGSRujarNOYWLqROg1ehyoOIDsymxVxiQi8hQGOUQeEBSgQ1JEoLyfVeLklJVo1hnbT9Wl5KK9w6jEUfL+0hxmc4ioa2OQQ+QhPWKUbE5WcY3X1OUIZ6adKY8Mcoioq2OQQ+QhPWIbg5wSF4IcN9TlOIKczYWbUVKnzjQYEZEnMMgh8pCezYqPnebY+Th3M2AxqXJeyaHJGBgzUO6V82v2r6qMSUTkCQxyiDw9XeVKJkc06wyOBawmJdBRydTuU+VxyaElqo1JRNTZGOQQeUhPx3RVca3zy8hFs87upyr3s9RrxzC1hxLkrM9fzykrIuqyGOQQebB/lYhRqk0WFFc7uYxc6HG6csxaod65haVhQPQATlkRUZfGIIfIQwINOiRHBLk+ZdXz9KPFxyrV5TTP5nDKioi6KgY5RF4wZeVS8XFcPyAkHrDUAYc3qHZu53Q/p2nKqriuWLVxiYg6C4McIi9o7+DSXjlizqvHaepPWYWnYWjsUDll9UPmD6qNS0TUWRjkEHX1FVbNp6wy1QtyhBm9Z8jjtwe+VXVcIqLOwCCHyIN6xSlBzsEilYKcw+sAcx3UMq3HNOi1euwu3Y19ZftUG5eIqDMwyCHyoPT4MHk8WFwDi9Xm2n45YcmAtQHIVq/FQ2RgJM5IPUPeX3xgsWrjEhH5XZCzfPlyzJgxA8nJydBoNPjqq688fUpEbpUSGYRAgxYNFhuyS51s1Omoy+mlBCM4+DvcMWW1+OBiWGwWVccmIvKbIKempgbDhg3DK6+84ulTIeoUWq0GfeJD5f2MwmrXBut9lnI8oG4rhtNTTkeUMQpFdUVYfni5qmMTEflNkDN9+nQ88cQTuOiiizx9KkSdPmW139Ugp5fSWBP524HqQqjFoDPgwvQL5f2FexeqNi4RkV8FOR1lMplQWVnZ4kbU1aQnKJmcfQVVrg0UGgckDVPuH/gNarqs72XQQIPVuauRVZGl6thERO7SpYOcp59+GhEREU23tLQ0T58SkdOZnIwCFzM5Qu+zleN+daesUsNSMTF1orzPbA4RdRVdOsiZO3cuKioqmm45OTmePiWiDktvrMk5UFQNq83JRp0Ofc4+msmxubBaqw1X9LtCHr/e/zVqzS4USRMRdZIuHeQYjUaEh4e3uBF1xUadRr0WJosNh8tcDB5SxwABoUBtMZC/DWqakDJBNu6sMldhUcYiVccmInKHLh3kEPkCnVaD3nGOuhwXp6z0AUc3Btz/M9Sk1Whx/aDr5f0FOxfAbDWrOj4RkU8HOdXV1diyZYu8CZmZmfJ+dna2p0+NqFOKjzMKq1QYTOkejr3q95u6oM8FiA2KRUFtAb7L/E718YmIfDbI2bBhA0aMGCFvwt133y3vP/zww54+NSK36pvQuIxcjeLjftOV45GNQGUe1GTUGXHtwGvl/Xd2vCObdxIReSuvCnImTZoEu93e6rZgwQJPnxqRWzk2BNynRiYnLBFIGa3c36d+NueyfpchLCAMmRWZ+DHzR9XHJyLyySCHyF/1T1QyOfvyq2F2pYdV04Dnum3KKsQQgusGXifvv7LlFdbmEJHXYpBD5AXSooIRZtSjwWpzfedjoV9jkHNwGWBSYbxjzBo4CzGBMcipysEXGV+oPj4RkRoY5BB5SQ+rgcnKFgg7jlS4PmBcfyCqJ2A1qd7LSgg2BOPWYbfK+29sfYP75hCRV2KQQ+QlBiVHyOPOXBXak4iu5P3PU+7v+hrucHHfi+W+OSX1JXhz25tu+QwiIlcwyCHyEoNTlEzOzlwVMjnCoIuO1uU01EBtBq0B955yr7z//s73kVGWofpnEBG5gkEOkZdlcnblVsLmansHIWWkMmUlWjC4oQBZmJQ2CWelnQWL3YIn/niCS8qJyKswyCHyEr3jQmR7h5oGK7JKatSZshp8sXJ/h/uKg+eOnYsgfRA2FW7CJ3s+cdvnEBF1FIMcIi+h12nRPylcvbocYcglyjHjZ6CuDO6QGJKIO0feKe8/t+E57Cnd45bPISLqKAY5RF5ksGOFlVp1OfEDgPhBgM0M7P4W7nJV/6swKXUSzDYz/rHsH1xtRURegUEOkRcZnHK0Lkc1jmzOlo/gLhqNBo9PeBzxwfHIqszCfSvug9VmddvnERG1B4McIi8yqDGTs/1IhWxpoophVwIaHZC9Bih031RSZGAknj/jeQRoA/B7zu/417p/qfdnICJyAoMcIi/SPzFcFh+X15pxoEilZd/hSUebdm50bx+44fHD8a/T/wUNNPhk7ydyo0AGOkTkKQxyiLxIgF6LYamR8v7GQ6XqDTzqBuW49WPAXAd3mtJ9StP+Oa9tfQ0vb36ZgQ4ReQSDHCIvM6pHlDxuyFJxNVTvM4GIbkB9udt2QG7umoHX4B+j/yHvz9s+D4/98RgbeRJRp9N3/kcS0YmM7q4EORsPqRjkaHXAyGuBpU8Aa98Ahl6u7KPjRtcOuhYBugA8tfYpfL7vcxwsP4jnzngOccFx8JQacw02F27GpoJNKKwtREVDBWAHAvWBMOqM8hgeEC7PMS4oTh5TQ1MRExTjsXMmIucxyCHyMqMag5yDxTUoqTYhJtSo0sDXAyueA3I3A1krgJ6nw92u6H8FkkOTcd/y++RmgRd+fSHuGX0PLuxzoVyR5W7VDdUyqFlfsB4b8jdgV8kuWO0dX/UlAp1h8cMwPG44RsSPQHpUOrQaJsKJvJ3G7kNVgZWVlYiIiEBFRQXCw5VVKkRd0eQXlmF/YTXemjUKUwclqjfwd38H1r8N9D4bmLUInSWzIlMGOrtLd8vHg2IG4c9D/yzbQqgVLIj/lOXV5GFnyU5sK9qmBDWlu1q1mkgJTcHohNHoEdEDEcYIWSRtsppQZ6lDvaUeFaYKFNUVKbfaIuTX5MMu0j3NRBmjMC5pHCZ3n4yJqRPljs9E5H3f3wxyiLzQ/V9swyfrc/CXM3ph7vQB6g1cmgm8PBIQX/x/WQEkDUVnsdgs+GDXB7IYWQQUjoDj3J7n4uxuZ6N/dH/oxLRaO8fKrszG3rK92Fu6V+6yLLI0ZabWU3y2hmhYanvBXtsDqXWhOCNSh/FxJgxLCEB8iA4wBAHBMUBkdyCuH2AMbfH+qoYqbC/ajq1FW7GlaAu2FG5BraW26eciwLko/SLMGjhL/nmIqOMY5HjwIhF1tk835ODez7fJ+pzPbxuv7uCf36j0shr0J+BS9y4pb0tJXQk+2P2B7HNVba5uej7MEIbBsYPlVFC3sG6yDkYEEGJ6qbKhEgU1BciuypZBzf7y/TL70opdC2t9Eqz1KdDWpmJYvQVn4RDGafeirz0TgRrzyU8wrr8ylSeyXaJgW99yulDs6iyCnt8P/44lWUtwpPqIfF6n0eHSvpfir8P/KvcMIqL2Y5DjwYtE1NkOFlXjrOeXIUCnxbZ/TkWgoX0ZjnbJ3w68MVFEBMAtS5Vu5R4gsjnLcpbhx6wfsS5vHarMVR16vwiA+kSmQ2tOwc6sEFSUJyDUFIbzjLtwdcR29K9eD731aMZFsBlCUBGUihxLJA5VadEAHUJgQrfAWvTRFyKgrqjlhxgjgIEzgNE3tXmdxBTZmrw1WLBjgTwKonBZBDqX9bsMei3LHonag0GOBy8SUWcTX57j//Ub8irq8d6NY3BGX5VXJC36M7BtoZKxuPYbt6+0OhnRAmJP2R7sKdmDfWX7ZB1MSX2JzNboNXoEG4JlywhRxNwvqh/6RffDwbxAPLF4N3KKyjFZuxHXBK7COPsWaJsXFoenAL3OBHpOBFJPAaJ6AlqlBiivog7vrT6EBaszUW9W6nZuGRWGv/UtRfDhFcCe74CqvKNjpYwGzrgPSJ/S5vUSgdq/1v8LGWUZ8nGfyD54cOyDGJ042u3Xj6irq2RNjucuEpEn3Pv5Vny64TBuPq0nHjp/oLqDlx0CXhkNWBuAq78A0iejqyiorMdji3dh//a1uEr3Ky7Ur0EEjk57IWEw0O9coP+5QNLwkwZwhZX1eH7JPizckCMfp0QG4dWrR2J4SrjSCmPTe8DOL5VrJaSNBaY/CyQPb7NWaFHGIrkBYrmpXD53Zf8rcdfIu2SgRkRtY5DTDgxyyJd8uzUXcz7ejH4JYfjpb25Y7v3Tg8CaV4DYfsCtK1rVnngbq82OD9ZkYdWSz3GN7WucrtveMmMjenSJW2wfp8ZffaAY93+xHdmltTDoNPi/8wdi1rjuylL36iJg9UvAurcBUTQtVoSN+TNw9sNAQEirscQKrRc3vogvMr6Qj3tF9MILk15A78jezl8AIh9WyUyO5y4SkSeU1TRg5BM/Q2zysPaBs5EQHqjuB9SVAa+cAtQUAWfcD5w5F97qUFEFFv3vZUwrX4gB2mz5nF2jhWbADGX/n55nKBseuqiy3ox7P9uGH3fmy8czhyXj6YuGIMTYWFtTmQcseVAp3BZi+gAXvw0kj2hzvNW5q/F/K/8PhXWFsobokVMfwXm9znP5PIl8TSWDHM9dJCJPueCVldh6uALPXToMl4xKVf8DdiwCPr8B0BqUbE68isvV1WBpwIZvX0f8llfRTVMgnzLrgqAbfR20424Donq4pR5q/spMPP3DHpk9So8PxfzrTkG3mGbTTft/Ab6eA1TlKtfurIeA8Xc01fscu5rsvhX3YW3eWvn48n6Xy95eYjdoInLv9ze37CTyYqc3Fhwv33fMqh+1iGXkfacDNjPw1e0yqPAK5jqY17yJ8mcGYfTWh2WAU6mNQMX4uTDcsxva6c+4JcARxPTUzRN74ZM/j0N8mBEZhdW48LVV2JDVrGFqn8nAbasAkUkS1+6XR4D3ZwJVSgaoObEU/s3Jb+IvQ/8iHy/cuxA3/HgDiuuK3XL+RHQUgxwiLzYxXQlyVu4vhs3mhs3JRb3Jec8DYl+X3E3KVIwnmeuBNa/C8uJQGH66F5HmQhTYI7G8198Qcu8uREy9HwhS2l642yk9ovHtnNMwJCUCpTUNuGreWny1WdkTRwqOBi77HzDzZUAUFYtWGW+eDhxSlpI3JzY5nD1iNl6f/LpcYr6teBuu+u4quZKMiNyHQQ6RFxvRLRJhgXr5JbspW8WGnc1FpAAXvaXcX/cWsP1zdDqbDdj2qVIj9NMD0NcW4rA9Fk9rbsH+K1fj9Gv/CV1gy52IO4Oog1r4l3GYOjABDVYb7lq4BS8s2SuntJqCRNH4VOweHT8QqC4A3jsf2Phem+OdlnIaPjz3Q3QP7y5bUFz7w7VYIZarE5FbMMgh8mIGnRZTBiTI+99vbz0Vopq+5wAT/67c//qvwKHV6DQHlwHzJgGLbgEqspFnj8Z95lvw94R3cP1dj2NCf8+2SggO0OONa0bJFhvCf3/bj4e+2tEysyZWdN30MzDoIsBmAb69A1j6lCjwaTWe6JklAh3RP0t0RZ/922x8vOfjzvwjEfkNBjlEXm76kCR5/HFH3tEMgjtMekCpz7HUAx9dDuRthVsV7AQ+uESpZcnbihoE4Vnz5TjT9DzCTr0RH/xlIpIivKPxpVarkT3EnvrTEJm8+XBtNu75fCss1mbNP0XPq0veAU7/h/J42TPAN7MBa+tWEqIx6FtT3sIFvS+QDUSfWvsU3tj6hnv/fon8EIMcIi83MT0WIQE65FbUy5VWbqPTA5e+C3SfAJgqgfcvALL/UP9zyrKARX8BXp8A7P8ZNo0e71mnYWL9i/gs+DK8deNEufmhyGJ5m6vGdsN/Lh8OnVaDRZuO4M6FW2BuHuiICEistDr/RWUvnc0fAB9fCZiabVbYyKAz4PEJj+P2YbfLx69ueRXPb3iegQ6RirzvvyJE1ILoW3VW45TVD9ubtRlwB9GR+8pPlBYGYh+d92Ye3RPGVdWFwPf/AF4eDWz7RPbOWmM8DWfWP4tHzNdi5IB0/HjnxKYVZd7qguEpePWqkXLDwO+25eG2Dzai3tyslYQw+kbgio8AfZAM5GS2qrbZ6qxmK7luG36bXFIuvLfrPTy65lHZ5oKIXMcgh6gLmD44UR6/d/eUlRAYDlz3DdDvPEB0+hZdy7+8DahT2hR0WOFuYPHdwEvDlMJmmxk5UeNwkeVJXFlxO4oNKXjyT4Mx79pRiAn17l2XHaYNTsRb146GUa/FL7sLccv7G1DXcExg0m86cP1iICgaOLIReGcaUNFsdVYzswbOwmPjH4NWo5W7JN+/4n7Z7ZyIXKOx+9AkMDcDJF9V22DByMd/lo0kv7htPEZ174Rl1CKb8NvjwMr/KB3LQ+KB8bOVLIUx7OS9sUSDy+2fKUvTG9XFD8eT9Zfhg0JljxuRtXnqT4ORGtU1+zqt3l+Mm9/fgNoGK8b0jMY715+CUMfuyA5Fe4H//QmoPAJEpAGzvgRi09sc76esn2SAI3pgnZl2Jp474zluGkh+oZI7HnvuIhF5g7s/3SLrQC4fnYZnLhnaeR+cvRb46jag9IDy2BAC9DoD6DERiOwGBEUqNSeiY3feFqWOp2jP0fdrdLD2m46v9Ofj/s3hEDM74YF62RtK7OIse0N1YRsPleL6d9ajymTB8LRIvHfDGEQEG1q+qDxHCXRKMoDgGOCaL47bCkIsKb9r6V1osDVgQsoE/GfSfxCoV7mlB5GXYZDjwYtE5A3WZZbisjfXIDhAh3UPTm6dMXAnsROyyMqsfFH5oj4ZjU5267YPuhDLDRPxyK8FyCqplT+aMjABT1w4WP1eXB60/XAFZr2zFuW1ZgxMCsf/bhrTeuqtphj44GIlEAwIBa78GOjZduPVP/L+wJxf56DeWo+xiWPx37P+yy7m5NMqmcnx3EUi8gZiZvnsF5bhYFGNbBp55Zhuntm0r2A7kPEzkLtZyd7UVyjTVyJDkTBYyVD0OgP7q/R49NtdWJGhtC+ICzPK7M2MoUldPnvTlj35lbjm7XUorjahR0ww3rtxDLrHHNOh3FQFfHIVkLkcEL2rLp4PDJzZ5ngbCzbi9l9uR62lFiPjR+LVs19FqAiOiHxQJYMcz10kIm/x1vIDeOr7PRiWGoGvZ58Gb60feumXDLy9MlM2uAzQaXHTxJ7465l9Ojf75AEHi6px7TvrcLisDjEhAbJGZ1haZOvWFYtuBnZ/qywzn/4sMOaWNsfbWrQVt/18G6rMVRgSO0S2hRB77BD5mko26CSii0amyqXLYr+cbYedXO3kRr/vLcTUF5fjzeUHZYAzeUAClvztdNw3rb/PBzhCr7hQLLp9PAanhKOkpgFXvPUHftujdE9vYggELn0PGHU9YLcB398D/PJom7sjD4sbhrfPeVsGNtuLt+PmJTejrN5N7T2IfBCXkBN1IbGhRpw/NFnef21pYyGwFxBTNHd8vBnXv7teZjGSIwIx/7rRePu60egRe8yUjY+LDwvEJ38+Va4cqzNbccv7G/HxuuyWL9LqgPP/A5zZ2BB15QtKcXcbuyMPjBmId855B9GB0dhTugc3/nQjO5gTtRODHKIu5rZJveXxx535yCio8vTpYPG2XEx5YRm+2ZoLrQa46bSe+PnuM3B24waG/khkrUSQJ1aPiYzW3EXbWzb2FERd0hn3AjNfUQq1t36stNMQdTvH6BvVF+9OexfxQfHYX74fN/x4AwpqjskQEVErDHKIupi+CWE4Z5ASQLz+u+eyOaIz+l8/2oTZH21GWa0Z/RPD8PVfT5PFxSF+MDV1MqItxb8vGYo7zurT1NjzH59va9kGQhg5S9ll2hAMHPgVWHAeUNU6gOkV0QsLpi1AUkgSsiqzcP2P1yO3Orez/jhEXRKDHKIuaPaZymZyX2/NRVZxTad//i+7CjD1xWWyrYHo43TH2en4ZvZpGJLKotjmxCqyu6f2k409xXX6fONh3LhgParqj5mW6jtV2R05OFZpjDp/ClC8v9V1TwtPk4FOamgqDlcfxnU/XofsymOmwoioCYMcoi5IBBNn9ouTUyGPLd7VaZ8rshBPLN4ld/ktrm5A34RQfHX7BNw9pS8C9PzPyYkae7597Wi5x5FYUn/pG2uQV1HX8kUpo4CblgBRPYHyQ0qgk7O+1VjJocky0OkR3gP5Nfkyo3Ow4qA7/rqJujz+V4moi1I6dWvw255CmVlxt9zyOlz+5hq5NFwQtTffzmH2pr3O7B+PhX8+Ve4XtCe/Cn96dTV25Va2fFFMb+Cmn5W9hupKgfdmAHt/aDVWQkiCrNHpE9kHRXVFskZnX9k+1/+SiXwMgxyiLqp3XChuOq2XvP/o4p2tO2GraOmeQpz73xXYlF2OsEA93pw1StbeGPU6t32mr2bgvrx9PPrEhyK/sh4Xvb4Kn27IaVmQHBoHXLcY6DMFsNQpmwduXNBqrNigWLnqakD0AJTWl8pVVztLdnbuH4jIyzHIIerC5pzVB0kRgcgprXPLtJXFasMzP+7BDQvWy5YFQ1Mj8P0dE3HOIKUrOnWcaEb6xa3jMTE9VjZcvffzbbj7062oMVmOvsjY2PZh+DXKXjrf3gksfarVXjpRgVGYN3UehsYORYWpArf8dAs2FRxtiErk7xjkEHVhYhXTs5cMlauRP1qbjUWbDqs2dk5pLa6c90fTCq7rx/fAZ7eeirTortkx3JuIBp6ikec/zuknC5K/3HwEM15Zid15zaavdAbggleA0+9VHi97BvhmjtJHrPlYxgi8OeVN2fpB7IwsNgxcfHBxJ/+JiLyTxt4iT9q1sa0D+av//LIP//klA4EGLT77y3iXVjmJ/yQsXJ+DxxfvQk2DVe758szFQ3He0CRVz5kU67NKMeejzXL6StRY3XFWOm6d1FsuQW+y4R3gu78rWZ20scqOyeEt/z7qLHV4YMUD+CX7F/n4tmG3yZsv9gkj31PJ3lWeu0hE3k6sshJTSsv3FSEiyIAPbx6LwSkdD3QKq+px/xfbZTGzMKZHNJ6/bBizN52w55CYtvplt1JALjqZiwxdi7/DfT8BX9wCmCqAkHjg0neBHi37l9nsNvxn03/w7o535eNze56LxyY8BqPumI7oRF6GQY4HLxJRVyD2XhFtFTYeKpOBzhvXjMKpvWPa9d4Giw3/++MQ/vtrBirqzLKppphKufG0nnI6hdxPZNDErtGPfLNT1j+J6/6X03thzlnpCApoLPAuOQAsnAUU7lR2SZ7yKHDqbGX35GY+3/c5nvzjSVjsFgyOGYwXz3wRiSGsoyLvxSDHgxeJqKuoNllw3TvrZKAjvvdmn9kHs8/qc9xVUCIwEhvUvbsqC9mltfK5QcnheOGy4eiXGNbJZ09CUZUJ//xmJ77bnicfp0YF4ZEZgzBlYGObjIZaYPFdwLaFyuP0c5TandD4Fhfwj7w/8Pff/47KhkpEGaPw7BnPYlzSOF5k8koMcjx4kYi6ktoGCx79ZhcWbshpaup51Zg0jOgehfgwo5waOVhUgxUZRVi1v0Q2kXS87u9T++LSUanQN68HIY/4cUc+Hvt2J3Ir6uXjs/vHy2CnW0ywsspq/dvATw8CVpOyU/IFrwL9prUY43DVYdz9+93YXbobWo0Wc0bMwU2Db2KdDnkdBjkevEhEXdG3W3PxxHe7UFBpOuHreseF4PoJPXHRiBT2nPLCgPXl3/bj7RUHYbba5a7SN4zvgdvP7COnJFGwC1h0C1CwQ3nD6BuBKY8rS9Ab1Vvq8eTaJ/HV/q/k4zNSz8Cj4x9FTFD7pjKJOgODHA9eJKKuSrRh+GlnPr7ZkoucsjoUVtYjOiRAFhKP7BYpd+EVRa5cgePd9hdWyymslfuL5ePIYIOs1blmXDcY7Wbgt8eBNa8oLw5PAc55Chh4QVOtjqj3+SLjCzy19imYbWZEB0bjsfGP4Yy0Mzz5xyJqwiCnHRjkEJGvEoHK0r2FePr7PcgorJbPpUUH4fZJfXDRyBQYDy0Hvr0DKG9s2Nn7LGD6v4FYpQu6sLd0L+5fcT/2lyvNPy/rexn+PvrvCBYd0Ik8iEGOBy8SEZG3ELtQi2LxF37eh8IqZSpS1FqJXmJXjYxD2IZXgJX/UWp1tAZgzC3AaX9rKkw2WU3476b/4v1d78vHoqP5Q+MewoSUCR79c5F/q/SXfXJee+01/Pvf/0ZeXh4GDRqE//znP5g4cWK73ssgh4j8qV7n43U5sl4nr7E4WfQVmzksGVf2MWPQ1qeg2f+z8mKRqRHBzvg7gZCYptVXD618CAW1yt48k7tNxpyRc9ArQumHRtSZ/CLIWbhwIWbNmiUDnQkTJuDNN9/E22+/jV27dqFbt24nfT+DHCLyN2KPo6+3HMEbyw7gQFFN0/PdooJwZ8/DmF40H8FFW5QnA0KBUdcDp9wERPdCjbkGr255FR/u/lBuJChWYM3sPVPulJwcmuy5PxT5nUp/CHLGjh2LkSNH4vXXX296bsCAAbjwwgvx9NNPn/T9DHKIyF/ZbHasOlAs+2CJ5ee1DY6u9HbMCNqOfwR8gW6mjMZnNNCkTwHG/BnofTYyKg7g5c0vY2nOUvlzvVaPi9MvxnUDr0NaeJoH/1TUpdntJ3zeDnvT725FZQViYmJ9N8hpaGhAcHAwPvvsM/zpT39qev7OO+/Eli1bsGzZslbvMZlM8tY8yElLS2NNDhHB36eyft5VgMXb8rDmQIncJFJ8pZyp3YLrdEswSbe16bU1gYmo6XM+wkZeiozQIPx3y8tYm7dW/kwDDc5MOxPXDrpWNgDlKjznVzmK/ajqxa3h6H1xlPcbrKi3WFHX7Gfy5w1WmBpMsJmqYW+og7WhDjZxs9QD5jrAYgIsddBYGhBgN0EnbqiDFvUI0JkQoGuAXmuCTivqsyywa6ywaaywaqywNB7NGhssGlvj0S7vN2jsMMsb0KABTFo7GjQaeV/cRNBggwb2xo22bfK3y3FTnldeo2j++qOva7w1rgC01lmx+7bdqn9/6+EliouLYbVakZDQuKtnI/E4Pz+/zfeI7M6jjz7aSWdIRNQ1BAfoccHwFHkThcpbD5djZUYJVh2Iwe1HRiPedBjX6H7BZbplCK/PR8iOt4EdbyMOcbgk7HT0TrwU64wHsL9mE37L+U3eBsYMxCV9L5H9sEIMIfAXItgQbTbK6xpQVmNGhTiKx43Pldc0HmvNLQIXEbCYGu9rbQ2IQSXiNOWI1VQgAjUI19QiDLUI09QiFDUI1NYiUFcDo64ONp0JZp0JDToLanR2lOu0qNBqUanVok6vgcmggSlYg3qNBiaNOGpRr9XArEozVi18iddkcnJzc5GSkoLVq1fj1FNPbXr+ySefxP/+9z/s2bOn1XuYySEi6ngz14NF1dh+pAK7sgsRnP07Bpb9iom29QjRHM2M19qN+FLXFwsjgpEVXgi7Vpn+0sKIBO1Y9Aw8C6mB/RESaEBIgA7BRr1yDNAjxHj0GBKgh1Gvla1FxGaG4ubufmiiTqmy3oyqegsq68zyfo1JBBwWGXyITJfIktQ2ZkvEY/Fz8brKeguq5Hss8rEY61iBMCEK1YjWVCFSU4VoiGP10aOmHEH6cuh1ldDoatGga0CZTocynRZlWp0MWkp1OpRrtcpzOp1KAcpRBmhghA4BGj0CNDoYNQYY5H1xM8CoNSBAG9DiZtQFymauBp0RRm0gAnRBMBoCYdQHI0AXiACDOAZBp9NDpxF5Pi20Gk3jDdBpdVD+FI6/Y/G8VmYAlVdroNEo7xHPaR33oUF1dQ2G9h/tu5mc2NhY6HS6VlmbwsLCVtkdB6PRKG9ERNQ+4ssnPSFM3i4amQpgJIC7UVZegdwt30GTsQTxhSsQbi7G1bbtuLoMKK3Q4qvQUHwaFoEjASbk2ZYjr3Y5DOUh0FX1R1XFGNSbxOKQ9n1R67UaGfgENAt+jj7WQq/VwmZXKjbE0SbnNZSjeFb801yModVq5FEQAY1oLisyKo5WJSdmhxFmhKIOIZp6hKNGBi3dG4OXKE01olCFKEOVfByqq4ZWVwOztg6VOhuK9VoU63TyVqLTIUunw2YRvGh1qGpqi2IA0KyT/EmIYCRSH4LIgDBEGSMQYYxCVFAsIoJjERUUg3BjBAJ1gQjUK8GIuIn7xz4nbiLg6EpEuYlPZ3IchcejRo2Sq6scBg4ciAsuuICFx0REnUV8LRTsBPb/DGvWGiB3E3S1RTLo2GIMwOdhofg5JBh12qNTG2lmG4bX69GrIQwRDTGos8Si0haAaose1VYdTDDAZDcoR3kLaHosC6GVDxaVG/K+ptl9LWzQwwq9xgqDOMKiPJbPW2DQiPviZxaEoF4GLSJ4idSZEKmrR5jWhFDxPGoRbK9DoLjZaqGFVU4B5ep1yNXrkafXo7Ax4yIyLWVa5Viq07b4s7aHOO8IfQiijJGICopGVHA8Io2RcrdpcYwKjGq6ycfGKATpg/y27qnSn5aQv/HGG3LK6q233sK8efOwc+dOdO/e/aTv5+oqIiI3EF8TVXlA7malX1blYdSVHcKK2hz8aK/C8kADTMcEATEWK/o1NKBfgxm9zWZ0N5vRw2xBpK319I871Gk0KNLpUKjXKUedDkX6xmOz+x0JXsR0j8iwxAXFITYkHrFBcbIHWGxQLGICY2QA4whcIgIiulw2xZPc9f3tNdNVwuWXX46SkhI89thjcjPAwYMH4/vvv29XgENERG4isgvhycqt/3nyqSAAU8XNbkdN5RGsOvgD1hZswPryPcgylaJEr8NqfRBWB4tXHhVp16K7uFmBWKsVQVYLgu12BNnFmBoENx01sq7EqgGsGg2sWj2sWh0sjbcqrRbVGg2qtBpUaoBijQ2FsKAIVhTZzaiyixVl7SOCk+SQZCSFJiEhOEE+dgQsjvviJgqu/TXT0lV5VSbHVczkEBF5Xq25VvbH2lu2V/bLyqrMwqHKQ8ivaXulrLuI6Z/44HiZeYkLjkN8ULxybHxOHoPj5OvIs/wik0NERF2faPg5NG6ovDVXZ6lDdmV2U9BTbiqXAZF4XtxqLc3um2tlx3S9Ri+nfXSaxptWJzcrDDOEITQgFGEBYfLmCGSaBzTMvBCDHCIi6hQiY9Ivup+8EXUG39r1h4iIiKgRgxwiIiLySQxyiIiIyCcxyCEiIiKf5FOFx47V8O7aHpqIiIjU5/jeVntXG58KcsRGgkJaWpqnT4WIiIic+B4X++WoxaeCnOjoaHnMzs5W9SL5a1QtgsWcnBxVN2byN7yOvJbeiL+XvI7eRmwC2K1bt6bvcbX4VJCjbexBIgIcfjGrQ1xHXkteR2/C30leS2/D30n1v8dVG0/V0YiIiIi8BIMcIiIi8kk+FeQYjUY88sgj8ki8lt6Av5O8lt6Iv5e8jv7yO+lTXciJiIiIfDKTQ0REROTAIIeIiIh8EoMcIiIi8kkMcoiIiMgndbkg57XXXkPPnj0RGBiIUaNGYcWKFSd8/bJly+TrxOt79eqFN954o9PO1Zeu5aJFizBlyhTExcXJja9OPfVU/PTTT516vr7yO+mwatUq6PV6DB8+3O3n6KvX0mQy4cEHH0T37t3lqozevXvjnXfe6bTz9aVr+eGHH2LYsGEIDg5GUlISbrjhhqZWOf5q+fLlmDFjBpKTk6HRaPDVV1+d9D38zlHnWqr2nWPvQj755BO7wWCwz5s3z75r1y77nXfeaQ8JCbEfOnSozdcfPHjQHhwcLF8nXi/eJ97/+eef2/1dR6+l+PkzzzxjX7dunX3fvn32uXPnyvdv2rTJ7s86eh0dysvL7b169bJPnTrVPmzYsE47X1+7ljNnzrSPHTvW/vPPP9szMzPta9euta9atcru7zp6LVesWGHXarX2l156Sf53UzweNGiQ/cILL7T7s++//97+4IMP2r/44guxCtn+5ZdfnvD1/M5R71qq9Z3TpYKcMWPG2G+99dYWz/Xv399+//33t/n6e++9V/68ub/85S/2cePG2f1dR69lWwYOHGh/9NFH7f7M2et4+eWX2x966CH7I488wiDHyWv5ww8/2CMiIuwlJSWu/SX6oI5ey3//+98y6G7uv//9rz01NdWt59mVtOeLmd856l1Ltb5zusx0VUNDAzZu3IipU6e2eF48Xr16dZvvWbNmTavXn3POOdiwYQPMZjP8lTPX8lg2mw1VVVWqN1Pzh+v47rvv4sCBA3LjK3L+Wn7zzTcYPXo0nn32WaSkpKBv37645557UFdX59eX1ZlrOX78eBw+fBjff/+9+IcvCgoK8Pnnn+O8887rpLP2DfzOcR9nv3O6TIPO4uJiWK1WJCQktHhePM7Pz2/zPeL5tl5vsVjkeGLe2R85cy2P9fzzz6OmpgaXXXYZ/JUz1zEjIwP333+/rI8Q9Tjk/LU8ePAgVq5cKWtOvvzySznG7bffjtLSUr+uy3HmWoogR9TkXH755aivr5f/jZw5cyZefvnlTjpr38DvHPdx9juny2RyHETBUnPiXx3HPney17f1vD/q6LV0+Pjjj/HPf/4TCxcuRHx8PPxde6+j+OK56qqr8Oijj8qsAzl/LR3/shM/E1/OY8aMwbnnnosXXngBCxYs8PtsTkev5a5du3DHHXfg4YcfllmgH3/8EZmZmbj11lv5a9pB/M5RnyvfOV3mn5KxsbHQ6XSt/iVSWFjY6l8sDomJiW2+XvwLOiYmBv7KmWvpIH7JbrrpJnz22WeYPHky/FlHr6NItYqp0s2bN2P27NlNX9Tiy0f8Ti5ZsgRnnXUW/JEzv5MiEyumqSIiIpqeGzBggLyeYuolPT0d/siZa/n0009jwoQJ+Mc//iEfDx06FCEhIZg4cSKeeOIJv816dxS/c9Tn6ndOl8nkBAQEyGWQP//8c4vnxWORam2LWHJ27OvFF4mYxzcYDPBXzlxLRzR9/fXX46OPPuJcvRPXUSyD3L59O7Zs2dJ0E/9S7tevn7w/duxY+CtnfifFl3Jubi6qq6ubntu3bx+0Wi1SU1Phr5y5lrW1tfK6NScCJYHtDduP3znqUuU7x94Fl0XOnz9fLou866675LLIrKws+XOxcmDWrFmtlvP97W9/k68X7+MScueu5UcffWTX6/X2V1991Z6Xl9d0E0uh/VlHr+OxuLrK+WtZVVUlV/9ccskl9p07d9qXLVtmT09Pt9988812f9fRa/nuu+/K/3+/9tpr9gMHDthXrlxpHz16tFyl5c/E79jmzZvlTXxdvvDCC/K+Yyk+v3Pcdy3V+s7pUkGOIP7A3bt3twcEBNhHjhwp/8PmcN1119nPOOOMFq///fff7SNGjJCv79Gjh/3111/3wFl3/Wsp7otfzGNv4nX+rqO/k80xyHHtWu7evds+efJke1BQkAx47r77bnttba2Kf7v+cy3FknGxRFdcy6SkJPvVV19tP3z4sN2fLV269IT/3eN3jvuupVrfORrxPypnmIiIiIg8rsvU5BARERF1BIMcIiIi8kkMcoiIiMgnMcghIiIin8Qgh4iIiHwSgxwiIiLySQxyiIiIyCcxyCEiIiKfxCCHiIiIfBKDHCIiIvJJDHKIyCuUlJQgPj4eWVlZbv2cSy65BC+88IJbP4OIvAN7VxGRV7jnnntQVlaG+fPnu/Vztm3bhjPPPBOZmZkIDw9362cRkWcxk0NEncpisbR6rq6uTgY3N998s9s/f+jQoejRowc+/PBDt38WEXkWgxwichsx9aTRaPD555/j9NNPh9FoxJdfftnqdT/88AP0ej1OPfXUpucmTZqEO+64A/feey+io6ORmJiIf/7zny3eJ14zZ84c3HXXXYiKikJCQgLeeust1NTU4IYbbkBYWBh69+4tx29u5syZ+Pjjj/k3T+TjGOQQkdts2bJFHp955hn83//9H3bu3ImpU6e2et3y5csxevToVs+/9957CAkJwdq1a/Hss8/isccew88//9zqNbGxsVi3bp0MeG677TZceumlGD9+PDZt2oRzzjkHs2bNQm1tbdN7xowZI19vMpnc8ucmIu/AIIeI3Gbr1q0ySPnss88wZcoU9OnTBxEREW1mfJKTk9ucWnrkkUeQnp6Oa6+9VgZCv/76a4vXDBs2DA899JB8zdy5cxEUFCSDnltuuUU+9/DDD8uiZlGL45CSkiIDnPz8fDf9yYnIGzDIISK3ZnLE1JCogTkRUZMTGBjYZpDTXFJSEgoLC4/7Gp1Oh5iYGAwZMqTpOTGFJTR/nwiEhObZHSLyPQxyiMitmRxRN3MyIvMiVlYdy2AwtHgs6ntsNttJX9P8OfFYaP6+0tJSeYyLi2v3n4WIuh4GOUTkFpWVlXIaasSIESd9rXjNrl27Ou1vYseOHUhNTZXBFRH5LgY5ROS2LI5Wq20xdXQ8ojhYFCW3lc1xhxUrVrRZAE1EvoVBDhG5Lcjp379/m7U2xxKBkCgq/vTTT93+t1FfXy+XsYvCZCLybdzxmIi8wvfffy93PRZTSSID5C6vvvoqvv76ayxZssRtn0FE3kHv6RMgIhLOPfdcZGRk4MiRI0hLS3PbRRFFyS+//DIvOpEfYCaHiIiIfBJrcoiIiMgnMcghIiIin8Qgh4iIiHwSgxwiIiLySQxyiIiIyCcxyCEiIiKfxCCHiIiIfBKDHCIiIvJJDHKIiIgIvuj/AQNp0GBkKrbDAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "r_plot = np.linspace(0, r_cut.value_in_unit(openmm.unit.nanometer), n_entries)\n",
    "plt.plot(r_plot, u_aa.value_in_unit(openmm.unit.kilojoule_per_mole), label=\"A-A\")\n",
    "plt.plot(r_plot, u_ab.value_in_unit(openmm.unit.kilojoule_per_mole), label=\"A-B\")\n",
    "plt.plot(r_plot, u_bb.value_in_unit(openmm.unit.kilojoule_per_mole), label=\"B-B\")\n",
    "plt.xlim(r_plot[0], r_plot[-1])\n",
    "plt.ylim(-0.2, 5.2)\n",
    "plt.xlabel(r\"$r\\ (\\rm{nm})$\")\n",
    "plt.ylabel(r\"$u(r)\\ (\\rm{kJ/mol})$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a78bf8cd-6faa-40ce-8c76-e7b59cf38d8a",
   "metadata": {},
   "source": [
    "## One Component\n",
    "\n",
    "Suppose for a moment that we wish to simulate a pure fluid of component A only.  Then it is straightforward to set up a [CustomNonbondedForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomNonbondedForce.html) to evaluate the potential:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3d268873-0795-4982-9f77-a4cfe1fd3c48",
   "metadata": {},
   "outputs": [],
   "source": [
    "function = openmm.Continuous1DFunction(u_aa, 0, r_cut)\n",
    "\n",
    "force = openmm.CustomNonbondedForce(\"table(r)\")\n",
    "force.addTabulatedFunction(\"table\", function)\n",
    "force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)\n",
    "force.setCutoffDistance(r_cut)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "331b417c-be7d-4129-854a-b048236db7f5",
   "metadata": {},
   "source": [
    "You can then add the force to a System, and call `force.addParticle()` for each particle in the System.\n",
    "\n",
    "There are a few things to keep in mind with this approach:\n",
    "\n",
    "* Internally, [Continuous1DFunction](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.Continuous1DFunction.html) uses natural cubic splines to interpolate between the table values.\n",
    "* Continuous1DFunction assumes that the table points are uniformly spaced between the lower and upper bounds.  If this is not the case for your tabulated values, you could resample them, or pass an appropriate function of `r` as the argument to `table` in the custom expression.\n",
    "* OpenMM converts tabulated function values to its [default unit system](https://docs.openmm.org/latest/userguide/theory/01_introduction.html#units), so results should be consistent as long as custom force expressions are dimensionally consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b73ae8f-c568-49d6-9d2e-c3c8196cccdc",
   "metadata": {},
   "source": [
    "## Multiple Components\n",
    "\n",
    "The simplest way to implement tabulated functions between particles of multiple types is to use a [Continuous3DFunction](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.Continuous3DFunction.html) in place of the Continuous1DFunction.  We can use the $x$ and $y$ dimensions of the function to look up potentials given a pair of particle types, and the $z$ dimension to look up potential values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2d120fee-9377-4cfc-8ee6-0c07cd054193",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_types = 2\n",
    "\n",
    "# For multidimensional tabulated functions, OpenMM expects values in column-major (\"Fortran\") order.\n",
    "# Ensure that the array is symmetric with respect to exchange of particle types.\n",
    "values = np.array([\n",
    "    [u_aa.value_in_unit(openmm.unit.kilojoule_per_mole), u_ab.value_in_unit(openmm.unit.kilojoule_per_mole)],\n",
    "    [u_ab.value_in_unit(openmm.unit.kilojoule_per_mole), u_bb.value_in_unit(openmm.unit.kilojoule_per_mole)],\n",
    "]).flatten(order=\"F\")\n",
    "\n",
    "function = openmm.Continuous3DFunction(n_types, n_types, n_entries, values, 0, n_types - 1, 0, n_types - 1, 0, r_cut)\n",
    "\n",
    "force = openmm.CustomNonbondedForce(\"table(type1, type2, r)\")\n",
    "force.addTabulatedFunction(\"table\", function)\n",
    "force.addPerParticleParameter(\"type\")\n",
    "force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)\n",
    "force.setCutoffDistance(r_cut)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "041d3ec4-fa33-404c-be29-c486e105e805",
   "metadata": {},
   "source": [
    "In this case, you must call `force.addParticle([particle_type])` to set the particle type index for each particle in the system.  (See the cookbook entry on [custom mixing rules](Customizing%20Lennard-Jones%20Mixing.ipynb) for another example of using particle types in OpenMM custom forces.)\n",
    "\n",
    "In general, Continuous3DFunction will perform interpolation along all of its dimensions.  However, in this example, it will only have an effect along the $z$ dimension as long as all particle types are integers, since the bounds of the function have been set so that the samples in the $x$ and $y$ dimensions fall on integer values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8cf5e3b-8b58-4adb-8b34-12cb1390abeb",
   "metadata": {},
   "source": [
    "## Notes\n",
    "\n",
    "Although this example shows a CustomNonbondedForce, you can also use tabulated functions with [CustomCentroidBondForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomCentroidBondForce.html), [CustomCompoundBondForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomCompoundBondForce.html), [CustomGBForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomGBForce.html), [CustomHbondForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomHbondForce.html), and [CustomManyParticleForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomManyParticleForce.html).  Note that [CustomBondForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomBondForce.html), [CustomAngleForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomAngleForce.html), and [CustomTorsionForce](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomTorsionForce.html) do not support tabulated functions, but their behavior can be fully replicated using CustomCompoundBondForce and the `distance(p1, p2)`, `angle(p1, p2, p3)`, and `dihedral(p1, p2, p3, p4)` functions supported in its custom energy expressions.\n",
    "\n",
    "Instead of 3D functions taking pairs of particle types along with distance, 2D functions taking a bonded interaction type and a distance or angle can be used.  Use [CustomCompoundBondForce.addPerBondParameter()](https://docs.openmm.org/latest/api-python/generated/openmm.openmm.CustomCompoundBondForce.html#openmm.openmm.CustomCompoundBondForce.addPerBondParameter) to add a type index to each bonded term, and use the type index to look up the potential values in tables.  Pay close attention to the boundary conditions associated with each kind of bonded term:\n",
    "* Typically, unlike pairwise potentials, bonds will not have cutoff distances beyond which the potential becomes 0.  You may wish to modify the custom energy expression using `select()` or another approach to define, *e.g.*, a linear extrapolation of the potential outside of the tabulated range.\n",
    "* Angles should be in $[0,\\pi]$.  Likewise, dihedrals should be in $[-\\pi,\\pi]$, and furthermore are periodic.  You can set `periodic=True` on OpenMM's continuous tabulated functions to ensure that there will be no discontinuity at the dihedral angle crossover point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f66e412-884d-455a-83bb-c32f31158f45",
   "metadata": {},
   "source": [
    "## Acknowledgements\n",
    "\n",
    "1. The tabulated potentials used in this tutorial were generated following [this OpenMSCG tutorial](http://software.rcc.uchicago.edu/mscg/tutorials/lesson-01/README.html) extended to a binary mixture.  Thanks to Clay Batton for providing these data."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.0"
  },
  "required_files": [
   "notebooks/cookbook/tabulated_potentials.dat"
  ]
 },
 "nbformat": 4,
 "nbformat_minor": 5
}